home *** CD-ROM | disk | FTP | other *** search
/ El Mac 9 / El Mac 9.iso / Shareware / Demos / Igor Demo Pro / 1 PutContentsIn Igor Pro Folder / WaveMetrics Procedures / Analysis / DSP (Fourier etc) / Power Spectral Density < prev   
Encoding:
Text File  |  1994-11-14  |  2.8 KB  |  91 lines  |  [TEXT/IGR0]
  1.  
  2. #include <DSP Window Functions>
  3. #include <BringDestToFront>
  4.  
  5. | Given a long data wave create a short result wave containing the Power
  6. | Spectral Density.  The name of the new wave is the name of the source + _psd
  7. | Note: if you don't want the baggage of all the window functions, create your own version
  8. | with your favorite and remove the window parameter (or choose a subset)
  9. | Also, note that a unity amplitude sine wave will give an integral peak size of
  10. | 0.5 rather than 0.707 as one would get if a real sine wave of 1 Volt peak amplitude
  11. | were applied across a 1 ohm resistor.  Thus if you want physical meaning for
  12. | the results you should scale by 1.414/R where R is the actual resistance.
  13. |
  14. |    Change 941114: Divided DC component by 2
  15. | 
  16. Macro PSD(w,seglen,window)
  17.     string w
  18.     Prompt w "data wave:",popup WaveList("*",";","")
  19.     variable seglen=1
  20.     Prompt seglen,"segment length:",popup "256;512;1024;2048;4096;8192"
  21.     variable window=2
  22.     Prompt window,"Window type:",popup "Square;Hann;Parzen;Welch;Hamming;"
  23.                         "BlackmanHarris3;KaiserBessel"
  24. ;
  25.     PauseUpdate; silent 1
  26.     
  27.     variable npsd= 2^(7+seglen)                | number of points in group (resultant psd wave len= npsd/2+1)
  28.     variable psdOffset= npsd/2                    | offset each group by this amount
  29.     variable psdFirst=0                            | start of current group
  30.     variable nsrc= numpnts($w)
  31.     variable nsegs,winNorm                        | count of number of segements and window normalization factor
  32.     string destw=w+"_psd",srctmp=w+"_tmp"
  33.     string winw=w+"_psdWin"                    | window goes here
  34.     
  35.     if( npsd > nsrc/2 )
  36.         Abort "psd: source wave should be MUCH longer than the segment length"
  37.     endif
  38.     make/o/n=(npsd/2+1) $destw
  39.     make/o/n=(npsd) $srctmp,$winw; $winw= 1
  40.     if( window==1 )
  41.         winNorm= 1
  42.     else
  43.         if( window==2 )
  44.             Hanning $winw;winNorm=0.372                |  winNorm is avg squared value
  45.         else
  46.             if( window==3 )
  47.                 winNorm= Parzen($winw)
  48.             else
  49.                 if( window==4 )
  50.                     winNorm= Welch($winw)
  51.                 else
  52.                     if( window==5 )
  53.                         winNorm= Hamming($winw)
  54.                     else
  55.                         if( window==6 )
  56.                             winNorm= BlackmanHarris3($winw)
  57.                         else
  58.                             if( window==7 )
  59.                                 winNorm= KaiserBessel($winw)
  60.                             else
  61.                                 Abort "unknown window index"
  62.                             endif
  63.                         endif
  64.                     endif
  65.                 endif
  66.             endif
  67.         endif
  68.     endif    | (kinda makes you wish we had elseif or switch constructs, huh?)
  69.  
  70.     Duplicate/O/R=[0,npsd-1] $w $srctmp; $srctmp *= $winw; fft $srctmp
  71.     CopyScales/P $srctmp, $destw
  72.     $destw= magsqr($srctmp)
  73.     psdFirst= psdOffset
  74.     nsegs=1
  75.     do
  76.         Duplicate/O/R=[psdFirst,psdFirst+npsd-1] $w $srctmp;   $srctmp *= $winw
  77.         fft $srctmp;   $destw += magsqr($srctmp);   psdFirst += psdOffset; nsegs+=1
  78.     while( psdFirst+npsd < nsrc )
  79.     winNorm= 2/(winNorm*nsegs*npsd^2);  $destw *= winNorm
  80.     $destw[0] /= 2
  81.  
  82.     KillWaves $srctmp,$winw
  83.     BringDestFront(destw)
  84.     if( numpnts($destw) <= 129 )
  85.         Modify mode($destw)=4,marker($destw)=19,msize($destw)=1
  86.     else
  87.         Modify mode($destw)=0
  88.     endif
  89. end
  90.  
  91.